home *** CD-ROM | disk | FTP | other *** search
/ MacAddict 84 / MacAddict_084_2003_08.iso / pc / Software / Audio & Music / Audacity 1.1.3.dmg / nyquist / dspprims.lsp < prev    next >
Lisp/Scheme  |  2002-09-16  |  13KB  |  457 lines
  1. ;; dspprims.lsp -- interface to dsp primitives
  2.  
  3. ;; ARESON - notch filter
  4. ;; 
  5. (defun areson (s c b &optional (n 0))
  6.   (multichan-expand #'nyq:areson s c b n))
  7.  
  8. (setf areson-implementations
  9.       (vector #'snd-areson #'snd-aresonvc #'snd-aresoncv #'snd-aresonvv))
  10.  
  11. ;; NYQ:ARESON - notch filter, single channel
  12. ;;
  13. (defun nyq:areson (signal center bandwidth normalize)
  14.   (select-implementation-1-2 areson-implementations 
  15.    signal center bandwidth normalize))
  16.  
  17.  
  18. ;; hp - highpass filter
  19. ;; 
  20. (defun hp (s c)
  21.   (multichan-expand #'nyq:hp s c))
  22.  
  23. (setf hp-implementations
  24.       (vector #'snd-atone #'snd-atonev))
  25.  
  26. ;; NYQ:hp - highpass filter, single channel
  27. ;;
  28. (defun nyq:hp (s c)
  29.   (select-implementation-1-1 hp-implementations s c))
  30.  
  31.  
  32. ;; comb-delay-from-hz -- compute the delay argument
  33. ;;
  34. (defun comb-delay-from-hz (hz caller)
  35.   (let (delay len)
  36.     ; convert from hz to delay, iterate over array if necessary
  37.     (cond ((arrayp hz)
  38.            (setf len (length hz))
  39.            (setf delay (make-array len))
  40.            (dotimes (i len)
  41.                (let ((h (aref hz i)))
  42.                  (cond ((numberp h)
  43.                         (setf (aref delay i) (/ (float h))))
  44.                        (t
  45.                         (error "bad argument type" h))))))
  46.           ((numberp hz)
  47.            (setf delay (/ (float hz))))
  48.           (t
  49.            (error (format nil "~A hz must be a number" caller) h)))
  50.     delay))
  51.  
  52.  
  53. ;; comb-feedback-from-decay -- compute the feedback argument
  54. ;;
  55. (defun comb-feedback (decay delay)
  56.   (let (feedback len d)
  57.     (cond ((arrayp decay)
  58.            (setf len (length decay))
  59.            (setf feedback (make-array len))
  60.            (dotimes (i len)
  61.                (setf d delay)
  62.                (cond ((arrayp d)
  63.                       (setf d (aref delay i))))
  64.                (setf (aref feedback i)
  65.                      (s-exp (scale (* -6.9078 d)
  66.                                    (recip (aref decay i)))))))
  67.           ((numberp decay)
  68.            (setf feedback (exp (/ (* -6.9078 delay) decay))))
  69.           (t
  70.            (setf feedback (s-exp (scale (* -6.9078 delay) (recip decay))))))
  71.     feedback))
  72.  
  73.  
  74. ;; COMB - comb filter
  75. ;; 
  76. ;; this is just a feedback-delay with different arguments
  77. ;;
  78. (defun comb (snd decay hz)
  79.   (let (delay feedback len d)
  80.     ; convert decay to feedback, iterate over array if necessary
  81.     (setf delay (comb-delay-from-hz hz "comb"))
  82.     (setf feedback (comb-feedback decay delay))
  83.     (feedback-delay snd delay feedback)))
  84.  
  85. ;; ALPASS - all-pass filter
  86. ;; 
  87. (defun alpass (snd decay hz)
  88.   (let (delay feedback len d)
  89.     ; convert decay to feedback, iterate over array if necessary
  90.     (setf delay (comb-delay-from-hz hz "alpass"))
  91.     (setf feedback (comb-feedback decay delay))
  92.     (nyq:alpass snd delay feedback)))
  93.  
  94.  
  95. ;; CONST -- a constant at control-srate
  96. ;;
  97. (defun const (value &optional (dur 1.0))
  98.   (let ((d (get-duration dur)))
  99.     (snd-const value *rslt* *CONTROL-SRATE* d)))
  100.  
  101.  
  102. ;; CONVOLVE - slow convolution
  103. ;; 
  104. (defun convolve (s r)
  105.   (multichan-expand #'snd-convolve s r))
  106.  
  107.  
  108. ;; FEEDBACK-DELAY -- (delay is quantized to sample period)
  109. ;;
  110. (defun feedback-delay (snd delay feedback)
  111.   (multichan-expand #'nyq:feedback-delay snd delay feedback))
  112.  
  113.  
  114. ;; NYQ:ALPASS -- alpass with delay and feedback args
  115. ;;
  116. (defun nyq:alpass (snd delay feedback)
  117.   (multichan-expand #'nyq:alpass1 snd delay feedback))
  118.  
  119.  
  120.  
  121. ;; SND-DELAY-ERROR -- report type error
  122. ;;
  123. (defun snd-delay-error (snd delay feedback)
  124.   (error "feedback-delay with variable delay is not implemented"))
  125.  
  126.  
  127. (setf feedback-delay-implementations
  128.       (vector #'snd-delay #'snd-delay-error #'snd-delaycv #'snd-delay-error))
  129.  
  130.  
  131. ;; NYQ:FEEDBACK-DELAY -- single channel delay
  132. ;;
  133. (defun nyq:feedback-delay (snd delay feedback)
  134.   (select-implementation-1-2 feedback-delay-implementations 
  135.                              snd delay feedback))
  136.  
  137.  
  138. ;; SND-ALPASS-ERROR -- report type error
  139. ;;
  140. (defun snd-alpass-error (snd delay feedback)
  141.   (error "allpass with variable hz is not implemented"))
  142.  
  143.  
  144. (if (not (fboundp 'snd-alpasscv))
  145.     (defun snd-alpasscv (snd delay feedback)
  146.       (error "snd-alpasscv (ALPASS with variable decay) is not implemented")))
  147.  
  148. (setf alpass-implementations
  149.       (vector #'snd-alpass #'snd-alpass-error 
  150.               #'snd-alpasscv #'snd-alpass-error))
  151.  
  152.  
  153.  
  154. ;; NYQ:ALPASS1 -- single channel alpass
  155. ;;
  156. (defun nyq:alpass1 (snd delay feedback)
  157.   (select-implementation-1-1 alpass-implementations 
  158.                              snd delay feedback))
  159.  
  160.  
  161. ;; S-EXP -- exponentiate a sound
  162. ;;
  163. (defun s-exp (s) (multichan-expand #'nyq:exp s))
  164.  
  165.  
  166. ;; NYQ:EXP -- exponentiate number or sound
  167. ;;
  168. (defun nyq:exp (s) (if (soundp s) (snd-exp s) (exp s)))
  169.  
  170.  
  171. ;; INTEGRATE -- integration
  172. ;;
  173. (defun integrate (s) (multichan-expand #'snd-integrate s))
  174.  
  175.  
  176. ;; S-LOG -- natural log of a sound
  177. ;;
  178. (defun s-log (s) (multichan-expand #'nyq:log s))
  179.  
  180.  
  181. ;; NYQ:LOG -- log of a number or sound
  182. ;;
  183. (defun nyq:log (s) (if (soundp s) (snd-log s) (log s)))
  184.  
  185.  
  186. ;; NOISE -- white noise
  187. ;;
  188. (defun noise (&optional (dur 1.0))
  189.   (let ((d (get-duration dur)))
  190.     (snd-white *rslt* *SOUND-SRATE* d)))
  191.  
  192.  
  193. (defun noise-gate (snd &optional (lookahead 0.5) (risetime 0.02) (falltime 0.5)
  194.                                                  (floor 0.01) (threshold 0.01))
  195.   (let ((rms (lp (mult snd snd) (/ *control-srate* 10.0))))
  196.     (setf threshold (* threshold threshold))
  197.     (mult snd (gate rms lookahead risetime falltime floor threshold))))
  198.  
  199.  
  200. ;; QUANTIZE -- quantize a sound
  201. ;;
  202. (defun quantize (s f) (multichan-expand #'snd-quantize s f))
  203.  
  204.  
  205. ;; RECIP -- reciprocal of a sound
  206. ;;
  207. (defun recip (s) (multichan-expand #'nyq:recip s))
  208.  
  209.  
  210. ;; NYQ:RECIP -- reciprocal of a number or sound
  211. ;;
  212. (defun nyq:recip (s) (if (soundp s) (snd-recip s) (/ (float s))))
  213.  
  214. ;; RMS -- compute the RMS of a sound
  215. ;;
  216. (defun rms (s &optional (rate 100.0) window-size)
  217.   (let (rslt step-size)
  218.     (setf step-size (round (/ (snd-srate s) rate)))
  219.     (cond ((null window-size)
  220.                (setf window-size step-size)))
  221.     (setf s (prod s s))
  222.     (setf result (snd-avg s window-size step-size OP-AVERAGE))
  223.         ;; compute square root of average
  224.         (s-exp (scale 0.5 (s-log result)))))
  225.  
  226.  
  227. ;; RESON - bandpass filter
  228. ;; 
  229. (defun reson (s c b &optional (n 0))
  230.   (multichan-expand #'nyq:reson s c b n))
  231.  
  232. (setf reson-implementations
  233.       (vector #'snd-reson #'snd-resonvc #'snd-resoncv #'snd-resonvv))
  234.  
  235. ;; NYQ:RESON - bandpass filter, single channel
  236. ;;
  237. (defun nyq:reson (signal center bandwidth normalize)
  238.   (select-implementation-1-2 reson-implementations 
  239.    signal center bandwidth normalize))
  240.  
  241.  
  242. ;; SHAPE -- waveshaper
  243. ;;
  244. (defun shape (snd shape origin)
  245.   (multichan-expand #'snd-shape snd shape origin))
  246.  
  247.  
  248. ;; SLOPE -- calculate the first derivative of a signal
  249. ;;
  250. (defun slope (s) (multichan-expand #'nyq:slope s))
  251.  
  252.  
  253. ;; NYQ:SLOPE -- first derivative of single channel
  254. ;;
  255. (defun nyq:slope (s)
  256.   (let* ((sr (snd-srate s))
  257.          (sr-inverse (/ sr)))
  258.     (snd-xform (snd-slope s) sr (- sr-inverse) 0.0 MAX-STOP-TIME 1.0)))
  259.  
  260.  
  261. ;; lp - lowpass filter
  262. ;; 
  263. (defun lp (s c)
  264.   (multichan-expand #'nyq:lp s c))
  265.  
  266. (setf lp-implementations
  267.       (vector #'snd-tone #'snd-tonev))
  268.  
  269. ;; NYQ:lp - lowpass filter, single channel
  270. ;;
  271. (defun nyq:lp (s c)
  272.   (select-implementation-1-1 lp-implementations s c))
  273.  
  274.  
  275.  
  276. ;;; fixed-parameter filters based on snd-biquad
  277.  
  278. (setf Pi 3.14159265358979)
  279.  
  280. (defun square (x) (* x x))
  281. (defun sinh (x) (* 0.5 (- (exp x) (exp (- x)))))
  282.  
  283.  
  284. ; remember that snd-biquad uses the opposite sign convention for a_i's 
  285. ; than Matlab does.
  286.  
  287. ; convenient biquad: normalize a0, and use zero initial conditions.
  288. (defun biquad (x b0 b1 b2 a0 a1 a2)
  289.   (let ((a0r (/ 1.0 a0)))
  290.     (snd-biquad x (* a0r b0) (* a0r b1) (* a0r b2) 
  291.                              (* a0r a1) (* a0r a2) 0 0)))
  292.  
  293. ; biquad with Matlab sign conventions for a_i's.
  294. (defun biquad-m (x b0 b1 b2 a0 a1 a2)
  295.   (biquad x b0 b1 b2 a0 (- a1) (- a2)))
  296.  
  297. ; two-pole lowpass
  298. (defun lowpass2 (x hz &optional (q 0.7071))
  299.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  300.          (cw (cos w))
  301.          (sw (sin w))
  302.          (alpha (* sw (sinh (/ 0.5 q))))
  303.          (a0 (+ 1.0 alpha))
  304.          (a1 (* -2.0 cw))
  305.          (a2 (- 1.0 alpha))
  306.          (b1 (- 1.0 cw))
  307.          (b0 (* 0.5 b1))
  308.          (b2 b0))
  309.     (biquad-m x b0 b1 b2 a0 a1 a2)))
  310.  
  311. ; two-pole highpass
  312. (defun highpass2 (x hz &optional (q 0.7071))
  313.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  314.          (cw (cos w))
  315.          (sw (sin w))
  316.          (alpha (* sw (sinh (/ 0.5 q))))
  317.          (a0 (+ 1.0 alpha))
  318.          (a1 (* -2.0 cw))
  319.          (a2 (- 1.0 alpha))
  320.          (b1 (- -1.0 cw))
  321.          (b0 (* -0.5 b1))
  322.          (b2 b0))
  323.     (biquad-m x b0 b1 b2 a0 a1 a2)))
  324.  
  325. ; two-pole bandpass.  max gain is unity.
  326. (defun bandpass2 (x hz q)
  327.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  328.          (cw (cos w))
  329.          (sw (sin w))
  330.          (alpha (* sw (sinh (/ 0.5 q))))
  331.          (a0 (+ 1.0 alpha))
  332.          (a1 (* -2.0 cw))
  333.          (a2 (- 1.0 alpha))
  334.          (b0 alpha)
  335.          (b1 0.0)
  336.          (b2 (- alpha)))
  337.     (biquad-m x b0 b1 b2 a0 a1 a2)))
  338.  
  339. ; two-pole notch.
  340. (defun notch2 (x hz q)
  341.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  342.          (cw (cos w))
  343.          (sw (sin w))
  344.          (alpha (* sw (sinh (/ 0.5 q))))
  345.          (a0 (+ 1.0 alpha))
  346.          (a1 (* -2.0 cw))
  347.          (a2 (- 1.0 alpha))
  348.          (b0 1.0)
  349.          (b1 (* -2.0 cw))
  350.          (b2 1.0))
  351.     (biquad-m x b0 b1 b2 a0 a1 a2)))
  352.  
  353. ; two-pole allpass.
  354. (defun allpass2 (x hz q)
  355.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  356.          (cw (cos w))
  357.          (sw (sin w))
  358.          (k (exp (* -0.5 w (/ 1.0 q))))
  359.          (a0 1.0)
  360.          (a1 (* -2.0 cw k))
  361.          (a2 (* k k))
  362.          (b0 a2)
  363.          (b1 a1)
  364.          (b2 1.0))
  365.     (biquad-m x b0 b1 b2 a0 a1 a2)))
  366.  
  367. ; bass shelving EQ.  gain in dB; Fc is halfway point.
  368. ; response becomes peaky at slope > 1.
  369. (defun eq-lowshelf (x hz gain &optional (slope 1.0))
  370.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  371.          (sw (sin w))
  372.          (cw (cos w))
  373.          (A (expt 10.0 (/ gain (* 2.0 20.0))))
  374.          (b (sqrt (- (/ (+ 1.0 (square A)) slope) (square (- A 1.0)))))
  375.          (apc (* cw (+ A 1.0)))
  376.          (amc (* cw (- A 1.0)))
  377.          (bs (* b sw))
  378.  
  379.          (b0 (*      A (+ A  1.0 (- amc)    bs  )))
  380.          (b1 (*  2.0 A (+ A -1.0 (- apc)        )))
  381.          (b2 (*      A (+ A  1.0 (- amc) (- bs) )))
  382.          (a0           (+ A  1.0    amc     bs  ))
  383.          (a1 (* -2.0   (+ A -1.0    apc         )))
  384.          (a2           (+ A  1.0    amc  (- bs) )))
  385.     (biquad-m x b0 b1 b2 a0 a1 a2)))
  386.  
  387. ; treble shelving EQ.  gain in dB; Fc is halfway point.
  388. ; response becomes peaky at slope > 1.
  389. (defun eq-highshelf (x hz gain &optional (slope 1.0))
  390.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  391.          (sw (sin w))
  392.          (cw (cos w))
  393.          (A (expt 10.0 (/ gain (* 2.0 20.0))))
  394.          (b (sqrt (- (/ (+ 1.0 (square A)) slope) (square (- A 1.0)))))
  395.          (apc (* cw (+ A 1.0)))
  396.          (amc (* cw (- A 1.0)))
  397.          (bs (* b sw))
  398.  
  399.          (b0 (*      A (+ A  1.0    amc     bs  )))
  400.          (b1 (* -2.0 A (+ A -1.0    apc         )))
  401.          (b2 (*      A (+ A  1.0    amc  (- bs) )))
  402.          (a0           (+ A  1.0 (- amc)    bs  ))
  403.          (a1 (*  2.0   (+ A -1.0 (- apc)        )))
  404.          (a2           (+ A  1.0 (- amc) (- bs) )))
  405.     (biquad-m x b0 b1 b2 a0 a1 a2)))
  406.  
  407. ; midrange EQ.  gain in dB, width in octaves (half-gain width).
  408. (defun eq-band (x hz gain width)
  409.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  410.          (sw (sin w))
  411.          (cw (cos w))
  412.          (J (sqrt (expt 10.0 (/ gain 20.0))))
  413.          (g (* sw (sinh (* 0.5 (log 2.0) width (/ w sw)))))
  414.          (b0 (+ 1.0 (* g J)))
  415.          (b1 (* -2.0 cw))
  416.          (b2 (- 1.0 (* g J)))
  417.          (a0 (+ 1.0 (/ g J)))
  418.          (a1 (- b1))
  419.          (a2 (- (/ g J) 1.0)))
  420.     (biquad x b0 b1 b2 a0 a1 a2)))
  421.  
  422. ; see failed attempt in eub-reject.lsp to do these with higher-order fns:
  423.  
  424. ; four-pole Butterworth lowpass
  425. (defun lowpass4 (x hz)
  426.   (lowpass2 (lowpass2 x hz 0.60492333) hz 1.33722126))
  427.  
  428. ; six-pole Butterworth lowpass
  429. (defun lowpass6 (x hz)
  430.   (lowpass2 (lowpass2 (lowpass2 x hz 0.58338080) 
  431.                                   hz 0.75932572) 
  432.                                   hz 1.95302407))
  433.  
  434. ; eight-pole Butterworth lowpass
  435. (defun lowpass8 (x hz)
  436.   (lowpass2 (lowpass2 (lowpass2 (lowpass2 x hz 0.57622191)
  437.                                             hz 0.66045510) 
  438.                                             hz 0.94276399)
  439.                                             hz 2.57900101))
  440.  
  441. ; four-pole Butterworth highpass
  442. (defun highpass4 (x hz)
  443.   (highpass2 (highpass2 x hz 0.60492333) hz 1.33722126))
  444.  
  445. ; six-pole Butterworth highpass
  446. (defun highpass6 (x hz)
  447.   (highpass2 (highpass2 (highpass2 x hz 0.58338080) 
  448.                                      hz 0.75932572) 
  449.                                      hz 1.95302407))
  450.  
  451. ; eight-pole Butterworth highpass
  452. (defun highpass8 (x hz)
  453.   (highpass2 (highpass2 (highpass2 (highpass2 x hz 0.57622191)
  454.                                                 hz 0.66045510) 
  455.                                                 hz 0.94276399)
  456.                                                 hz 2.57900101))
  457.